function result = sdm_panel_rdd_g(y,x,W,r,D,c,T,ndraw,nomit,prior)
% PURPOSE: MCMC SDM Regression discontinuity design model estimates 
%          for static spatial panels 
%          (N regions*T time periods) with spatial fixed effects (sfe) 
%          and/or time period fixed effects (tfe)
%          y = rho*W*y + X*b + D*tr + sfe(optional) + tfe(optional) + e, 
%          e = N(0,sige*I_NT), 
%          b   = N(c,C), default c = 0, C = eye(k)*1e+12
%          sige = gamma(nu,d0), default nu=0, d0=0   
%          no prior for rho
% Supply data sorted first by time and then by spatial units, so first region 1,
% region 2, et cetera, in the first year, then region 1, region 2, et
% cetera in the second year, and so on
% sar_panel_rdd_g with prior.model=1,2,3
%                transforms y and x to deviation of the spatial and/or time means
% sar_panel_rdd_g with prior.model=0
%                no transformation (recommended)
% ---------------------------------------------------
%  USAGE: results = sdm_panel_rdd_g(y,x,W,r,c,N,T,ndraw,nomit,prior)
%  where:  y = N*T x 1 dependent variable vector
%          x = N*T x k independent variables matrix
%          W = spatial weights matrix (standardized)
%          r = continuous or discrete treatment vector (N*T x 1)
%          c = scalar cut-off/treatment break, e.g., 0
%          N. B. W-matrix can be N*T x N*T or N x N
%          T = number of points in time
%       prior = a structure variable with input options:
%       prior.plot_flag = 1, plotting MCMC sampling
%                       = 0, no plotting (default)
%       prior.model = 0 pooled model without fixed effects (default, x may contain an intercept)
%                  = 1 spatial fixed effects (x may not contain an intercept)
%                  = 2 time period fixed effects (x may not contain an intercept)
%                  = 3 spatial and time period fixed effects (x may not contain an intercept)
%       prior.fe    = report fixed effects and their t-values in prt_panel
%                     (default=0=not reported; prior.fe=1=report) 
%               prior.beta, prior means for beta,   b (default (k x 1) vector = 0)
%               priov.bcov, prior beta covariance, T above (default eye(k)*1e+12)
%               prior.nu,   informative Gamma(nu,d0) prior on sige
%               prior.d0    informative Gamma(nu,d0) prior on sige
%                           default for above: nu=0,d0=0 (diffuse prior)
%       prior.rmin  = (optional) minimum value of rho to use in search  
%       prior.rmax  = (optional) maximum value of rho to use in search    
%       prior.lflag = 0 for full lndet computation (default = 1, fastest)
%                  = 1 for MC lndet approximation (fast for very large problems)
%                  = 2 for Spline lndet approximation (medium speed)
%       prior.order = order to use with info.lflag = 1 option (default = 50)
%       prior.iter  = iterations to use with info.lflag = 1 option (default = 30)  
%       prior.lndet = a matrix returned in results.lndet containing log-determinant information to save time
% ---------------------------------------------------
%  RETURNS: a structure
%         results.meth  = 'sdmrdd_g if prior.model = 0
%                       = 'sdmrdd_sfe_g' if prior.model=1
%                       = 'sdmrdd_tfe_g' if prior.model=2
%                       = 'sdmrdd_stfe_g' if prior.model=3
%         results.beta  = bhat (includes theta)
%         results.rho   = rho 
%         results.tr    = mean of treatment effects
%         results.lower = beta from lower kernel regression
%         results.upper = beta from upper kernel regression
%         results.bdraw = (ndraw-nomit)xk matrix of MCMC draws for beta
%         results.pdraw = (ndraw-nomit)x1 vector of MCMC draws for rho
%         results.trdraw = N*T x 1 vector of MCMC draws for tr (treatment effect)
%         results.samples= ndraw x 2 matrix of lower and upper bandwidth sample sizes
%         results.sdraw = (ndraw-nomit)x1 vector of MCMC draws for sige
%         results.bmean = b prior means (prior.beta from input)
%         results.bstd  = b prior std deviation, sqrt(diag(prior.bcov))
%         results.nu    = prior nu-value for sige prior (default = 0)
%         results.d0    = prior d0-value for sige prior (default = 0)
%         results.iprior = 1 for informative prior on beta, 
%                        = 0 for default no prior on beta
%         results.direct   = nvar x 5 matrix with direct effect, t-stat, t-prob, lower01, upper99
%         results.indirect = nvar x 5 matrix with indirect effect, t-stat, t-prob, lower01, upper99
%         results.total    = nvar x 5 matrix with total effect, t-stat, t-prob, lower01, upper99
%         results.direct_treat   = nvar x 5 matrix with direct treatment effect, t-stat, t-prob, lower01, upper99
%         results.indirect_treat = nvar x 5 matrix with indirect treatment effect, t-stat, t-prob, lower01, upper99
%         results.total_treat    = nvar x 5 matrix with total treatment effect, t-stat, t-prob, lower01, upper99
%         results.direct_draws   = ndraw x nvar matrix of direct effect draws
%         results.indirect_draws = ndraw x nvar matrix of indirect effect draws
%         results.total_draws    = ndraw x nvar matrix of total effect draws
%         results.direct_treat_draws   = ndraw x 1 vector of direct treatment effect draws
%         results.indirect_treat_draws = ndraw x 1 vector of indirect treatment effect draws
%         results.total_treat_draws    = ndraw x1 vector of total treatment effect draws
%         results.cov   = asymptotic variance-covariance matrix of the parameters b(eta) and rho
%         results.resid = y-p*W*y-x*b
%         results.sige  = (y-p*W*y-x*b)'*(y-p*W*y-x*b)/n
%         results.rsqr  = rsquared
%         results.corr2 = goodness-of-fit between actual and fitted values
%         results.sfe   = spatial fixed effects (if prior.model=1 or 3)
%         results.tfe   = time period fixed effects (if prior.model=2 or 3)
%         results.tsfe  = t-values spatial fixed effects (if prior.model=1 or 3)
%         results.ttfe  = t-values time period fixed effects (if prior.model=2 or 3)
%         results.con   = intercept 
%         results.con   = t-value intercept
%         results.nobs  = # of observations
%         results.nvar  = # of explanatory variables in x 
%         results.tnvar = # fixed effects
%         results.iter  = # of iterations taken
%         results.rmax  = 1/max eigenvalue of W (or rmax if input)
%         results.rmin  = 1/min eigenvalue of W (or rmin if input)
%         results.lflag = lflag from input
%         results.fe    = fe from input
%         results.liter = info.iter option from input
%         results.order = info.order option from input
%         results.limit = matrix of [rho lower95,logdet approx, upper95] intervals
%                         for the case of lflag = 1
%         results.time1 = time for log determinant calcluation
%         results.time2 = time for eigenvalue calculation
%         results.time4 = time for MCMC sampling
%         results.time  = total time taken      
%         results.lndet = a matrix containing log-determinant information
%                          (for use in later function calls to save time)
% --------------------------------------------------
% Written by Gary Cornwall BEA,
% modified by J.P. LeSage

%Inputs
    %y: an Nx1 vector of outcomes
    %x: an NxK matrix of information pertinent to y
    %w: an NxN adjacency matrix
    %r: an Nx1 vector of values for the running variable
    %c: a scalar threshold value that determines treatment
    %ndraw: number of iterations through the sampler
    %nomit: number of burn-in iterations
    %prior: each prior is for a different stage of the sampler, default is
            %an uniformative -- but proper -- prior.
            
    %outputs:
    
    
%set up covariates
% if sum(x(:,1)) == size(x,1), x = x(:,2:end); end
% pp = size(x,2); 
% x = [ones(size(x,1),1) x];
%Checking sizes
[Ny, ~] = size(y);
[Nx, k] = size(x);
x = [x W*x];

if Ny == Nx; else; disp('Inputs are not the rights size'); end

fe=0;
model=0;

k = size(x,2);

results.nvar = k;

fields = fieldnames(prior);
nf = length(fields);
novi_flag = 0;
results.rval = 0;
nu = 0; d0 = 0; % default to a diffuse prior on sige
cc = zeros(k,1); 
Tj = eye(k)*1e+12;
Q = inv(Tj);
Qpc = Q*cc;
iprior = 0;
results.iprior = 0;
for i=1:nf
    if  strcmp(fields{i},'nu')
        nu = prior.nu;
        results.nu = nu;
    elseif strcmp(fields{i},'d0')
        d0 = prior.d0;
        results.d0 = d0;
    elseif strcmp(fields{i},'plt_flag')
        plt_flag = prior.plt_flag;
    elseif strcmp(fields(i),'bw')
        if prior.bw == 1
            bw_calc ='IK';
        elseif prior.bw == 2
            bw_calc = 'CCT';
        end
    elseif strcmp(fields{i},'beta')
        c = prior.beta;
        if size(c,1) ~= k
            error('sarrdd_panel_FE_g: wrong size prior means, must be k x 1 vector');
        end
        results.iprior = 1;
        results.bmean = c;
    elseif strcmp(fields{i},'bcov')
        TI = prior.bcov;
        if size(TI,2) ~= k
            error('sarrdd_panel_FE_g: wrong size prior variance-covariance, must be k x k matrix');
        end
        results.iprior = 1;
        results.bstd = diag(sqrt(TI));
        
        Q = inv(TI);
        Qpc = Q*cc;
        
        
    end
end


fields = fieldnames(prior);
nf = length(fields);
if nf > 0
    for i=1:nf
        if strcmp(fields{i},'model') model = prior.model;
        elseif strcmp(fields{i},'fe') fe = prior.fe;
        end
    end
end
if model==0
    results.meth='sarrdd_g';
elseif model==1
    results.meth='sarrdd_sfe_g';
elseif model==2
    results.meth='sarrdd_tfe_g';
elseif model==3
    results.meth='sarrdd_stfe_g';
else
    error('sarrdd_panel_FE_g: wrong input number of prior.model');
end



% check size of user inputs for comformability
[nobs, nvar] = size(x);
N = nobs/T;

    [n1,n2] = size(W);
    
    if n1 == N*T
        Wlarge = 1;
    else
        Wlarge = 0;
        W = kron(eye(T),W);
    end

[nchk junk] = size(y);
if nchk ~= nobs
error('sarrdd_panel_FE_g: wrong size vector y or matrix x');
end;

% demeaning of the y and x variables, depending on (info.)model

Wy = W*y;

if (model==1 | model==3);
meanny=zeros(N,1);
meannwy=zeros(N,1);
meannx=zeros(N,nvar);
for i=1:N
    ym=zeros(T,1);
    wym=zeros(T,1);
    xm=zeros(T,nvar);
    for t=1:T
        ym(t)=y(i+(t-1)*N,1);
        wym(t)=Wy(i+(t-1)*N,1);
        xm(t,:)=x(i+(t-1)*N,:);
    end
    meanny(i)=mean(ym);
    meannwy(i)=mean(wym);
    meannx(i,:)=mean(xm);
end
clear ym wym xm;
end % if statement

if ( model==2 | model==3)
meanty=zeros(T,1);
meantwy=zeros(T,1);
meantx=zeros(T,nvar);
for i=1:T
    t1=1+(i-1)*N;t2=i*N;
    ym=y([t1:t2],1);
    wym=Wy([t1:t2],1);
    xm=x([t1:t2],:);
    meanty(i)=mean(ym);
    meantwy(i)=mean(wym);
    meantx(i,:)=mean(xm);
end
clear ym wym xm;
end % if statement
    
en=ones(T,1);
et=ones(N,1);
ent=ones(nobs,1);

if model==1
    ywith=y-kron(en,meanny);
    wywith=Wy-kron(en,meannwy);
    xwith=x-kron(en,meannx);
elseif model==2
    ywith=y-kron(meanty,et);
    wywith=Wy-kron(meantwy,et);
    xwith=x-kron(meantx,et);
elseif model==3
    ywith=y-kron(en,meanny)-kron(meanty,et)+kron(ent,mean(y));
    wywith=Wy-kron(en,meannwy)-kron(meantwy,et)+kron(ent,mean(Wy));
    xwith=x-kron(en,meannx)-kron(meantx,et)+kron(ent,mean(x));
else
    ywith=y;
    wywith=W*y;
    xwith=x;
end % if statement

y = ywith;
x = xwith;
wy = wywith;

%Setting up priors and initial conditions
%rho
p1 = 0.1; p2 = 0.1;
% p3 = 0.1;
%Variance
s1 = 1; s2p = 1; s2m = 1; 
s2 = 1;

% s3 = 1;
pnu2 = 1;
mnu2 = 1;
%set domain of rho
rmin = -0.999; rmax = 0.999;
[detval, ~] = sar_lndet(1,W,rmin,rmax,0,50,30);
%Stage 1 priors
%     if exist('prior1')==1, else, prior1 = struct(); end
%     if isfield(prior1,'bcov')==1, iB01 = prior1.bcov; else, iB01 = inv(diag(ones(k,1)*1e+12)); end
%     if isfield(prior1,'bmean')==1, b01 = prior1.bmean; else, b01 = zeros(k,1); end
%     if isfield(prior1,'a1')==1, a11 = prior1.a1; else, a11 = 1.01; end
%     if isfield(prior1,'a2')==1, a21 = prior1.a2; else, a21 = 1.01; end
%     if isfield(prior1,'nu')==1, nu1 = prior1.nu; else, nu1 = 0; end
%     if isfield(prior1,'d0')==1, d01 = prior1.d01; else, d01 = 0; end
%     TI = inv(iB01);
%     TIc = TI*b01;
% 
% %Stage 2 priors
%     if exist('prior2')==1, else, prior2 = struct(); end
%     if isfield(prior2,'mbcov')==1,  miB02 = prior2.mbcov; else, miB02 = inv(diag(ones(2,1)*1e+12)); end
%     if isfield(prior2,'pbcov')==1,  piB02 = prior2.pbcov; else, piB02 = inv(diag(ones(2,1)*1e+12)); end
%     if isfield(prior2,'mbmean')==1, mb02 = prior2.mbmean; else, mb02 = zeros(2,1); end
%     if isfield(prior2,'pbmean')==1, pb02 = prior2.pbmean; else, pb02 = zeros(2,1); end
%     if isfield(prior2,'md0')==1,    md02 = prior2.md0; else, md02 = 0; end
%     if isfield(prior2,'pd0')==1,    pd02 = prior2.pd0; else, pd02 = 0; end
%     if isfield(prior2,'mnu')==1,    mnu2 = prior2.mnu; else, mnu2 = 0; end
%     if isfield(prior2,'pnu')==1,    pnu2 = prior2.pnu; else, pnu2 = 0; end
%        % The distribution then has mean A*B and variance A*B^2.
% %Stage 3 priors        
%     if exist('prior3')==1, else, prior3 = struct(); end
%     if isfield(prior3,'bcov')==1, iB03 = prior3.bcov; else, iB03 = inv(diag(ones(k,1)*1e+12)); end
%     if isfield(prior3,'bmean')==1, b03 = prior3.bmean; else, b03 = zeros(k,1); end
%     if isfield(prior3,'a1')==1, a13 = prior3.a1; else, a13 = 1.01; end
%     if isfield(prior3,'a2')==1, a23 = prior3.a2; else, a23 = 1.01; end
%     if isfield(prior3,'nu')==1, nu3 = prior3.nu; else, nu3 = 1; end
%     if isfield(prior3,'d0')==1, d03 = prior3.d01; else, d03 = 10; end
%Storage Vectors
    beta1  = zeros(ndraw,k);
    sige1  = zeros(ndraw,1);
    rho1   = zeros(ndraw,1);
    mbeta2 = zeros(ndraw,2);
    pbeta2 = zeros(ndraw,2);
    msige2 = zeros(ndraw,1);
    psige2 = zeros(ndraw,1);
    bw_left  = zeros(ndraw,1);
    bw_right = zeros(ndraw,1);
%     beta3  = zeros(ndraw,k+1);
%     sige3  = zeros(ndraw,1);
    treat_save = zeros(ndraw,1);

    samples = zeros(ndraw,2);
    
    wxp_save = zeros(Nx,1);
    wxm_save = zeros(Nx,1);

% This is just used on the first trip through MCMC
%     wy = W*y;
    
    xpx = x'*x;
    xpxi = inv(xpx);
    xpy = x'*y;
  
    noo = 100;

%Sampler
% hwait = waitbar(0,'SDMRD: MCMC Sampling ...');
for j = 1:ndraw
    %Stage 1 Beta
    
          ys1 = y - p1*wy;          
          b0 = xpxi*x'*ys1;
          b1 = norm_rnd(s1*xpxi) + b0;  
          
%     ys1 = y - p1*wy;
% %     Db1 = (xpx.*s1 + iB01);
%      Db1 = xpx.*s1;
%     
% %     db1 = x'*ys1.*s1 + iB01*b01;
% 
%     H = chol(inv(Db1));
%     b1 = Db1\db1 + H'*randn(k,1);
    %Stage 1 Sige
    nu_1 = nu + Nx;
    d1_1 = d0 + (ys1 - x*b1)'*(ys1 - x*b1);
    chi_1 = chis_rnd(1,nu_1);
    s1 = d1_1/chi_1;
    %update rho
          b0 = xpxi*xpy;
          bd = xpxi*x'*wy;
          e0 = y - x*b0;
          ed = wy - x*bd;
          epe0 = e0'*e0;
          eped = ed'*ed;
          epe0d = ed'*e0;

%     b0_1 = (xpx + s1*iB01)\(xpy + s1*(iB01*b01));
%     bd_1 = (xpx + s1*iB01)\(x'*wy + s1*(iB01*b01));
%     e0_1 = y - x*b0_1;
%     ed_1 = wy - x*bd_1;
%     epe0_1 = e0_1'*e0_1;
%     eped_1 = ed_1'*ed_1;
%     epe0d_1 = ed_1'*e0_1;
    p1 = draw_rho(detval, epe0,eped,epe0d,Nx,k,p1);
    %Filter Residual
    F = speye(Nx) - p1*W;
    e_use = y - F\(x*b1);
%     size(e_use)
    Fe = F*e_use;
%     size(Fe)
    %Bandwidth
%     switch bw_calc
%         case 'IK'

%              [bw,fl] = ikbw(Fe,r,c);
% %             bw = 14.4509;
%             bw_l = bw; bw_r = bw;
%         case 'CCT'
%             [bw,bias] = rdbwselect_CCT(Fe,r,c);
%             bw_l = bw; bw_r = bw;
%         case 'AI'
%             AI_calc = rddmmse(r,Fe,c);
%             bw_l = AI_calc.h(1); bw_r = AI_calc.h(2);
%         otherwise
%             disp('BW calculation method not chosen, default to IK')
%             bw = ik_optbw(Fe,r,c);
%             bw_l = bw; bw_r = bw;
%     end

    %Bandwidth
    switch bw_calc
        case 'IK'
            [bw,flag] = ikbw(Fe,r,c);
            bw_l = bw; bw_r = bw;
        case 'CCT'
            [bw,bias] = rdbwselect_CCT(Fe,r,c);
            bw_l = bw; bw_r = bw;
        case 'AI'
            AI_calc = rddmmse(r,Fe,c);
            bw_l = AI_calc.h(1); bw_r = AI_calc.h(2);
        otherwise
            disp('BW calculation method not chosen, default to IK')
            bw = ik_optbw(Fe,r,c);
            bw_l = bw; bw_r = bw;
    end
    %Weights for observations
    [wyp, wym, wxp, wxm, ~, ~] = weights(Fe,r,c,bw_l,bw_r);
    
    wxp_save = wxp;
    wxm_save = wxm;
    %Regression from above cutoff
    xpx2 = wxp'*wxp;
    xpx2i = xpx2\eye(size(xpx2));
    xpy2 = wxp'*wyp;
    pbeta = xpx2i*xpy2;
    pbeta = norm_rnd(s2p*xpx2i) + pbeta;
%     pDb = wxp'*wxp.*s2p + piB02;
%     pdb   = wxp'*wyp.*s2p + piB02*pb02;
%     pH    = chol(inv(pDb));
%     pbeta = pDb\pdb + pH'*randn(2,1);
    %variance from above cutoff
    if  length(wyp) == 0
        pnu2 = samples(j-1,1);
    end
    lpnu  = pnu2 + length(wyp);
    
    pe    = wyp - wxp*pbeta;
    pd1   = pe'*pe;
    chi   = chis_rnd(1,lpnu);
    s2p   = pd1/chi;
    %Regression from below cutoff
    xpx3 = wxm'*wxm;
    xpx3i = xpx3\eye(size(xpx3));
    xpy3 = wxm'*wym;
    mbeta = xpx3i*xpy3;
    mbeta = norm_rnd(s2m*xpx3i) + mbeta;

%     mDb   = wxm'*wxm.*s2m + miB02;
%     mdb   = wxm'*wym.*s2m + miB02*mb02;
%     mH    = chol(inv(mDb));
%     mbeta = mDb\mdb + mH'*randn(2,1);    
    %Variance from below cutoff
    if  length(wym) == 0
        mnu2 = samples(j-1,2);
    end

    lmnu = mnu2 + length(wym);
    
    samples(j,:) = [length(wyp) length(wym)];
  
%     plot(samples);
%     drawnow;

    me   = wym - wxm*mbeta;
    md1  = me'*me;
    chi  = chis_rnd(1,lmnu);
    s2m  = md1/chi;
    
    treat = (pbeta(1,1) - mbeta(1,1));
    
        
% =============================================
    %Stage 3

    % update x to be [x D*treat]

    % update beta
    x2 = [xwith D*treat];
    xpx2 = x2'*x2;
    xpxi2 = inv(xpx2);
    xpy2 = x2'*y;

          ys1 = y - p2*wy;          
          b0 = xpxi2*x2'*ys1;
          b2 = norm_rnd(s2*xpxi2) + b0;  
          
    %Stage 3 
    % update Sige
    nu_1 = nu + Nx;
    d1_1 = d0 + (ys1 - x2*b2)'*(ys1 - x2*b2);
    chi_1 = chis_rnd(1,nu_1);
    s2 = d1_1/chi_1;
    %update rho
          b0 = xpxi2*xpy2;
          bd = xpxi2*x2'*wy;
          e0 = y - x2*b0;
          ed = wy - x2*bd;
          epe0 = e0'*e0;
          eped = ed'*ed;
          epe0d = ed'*e0;
    p2 = draw_rho(detval, epe0,eped,epe0d,Nx,k,p2);
    %Update beta
%     ys3 = y-p3*wy - treat*double(r>c);
%     Db3 = xpx.*s3 + iB03;
%     db3 = x'*ys3.*s3 + iB03*b03;
%     H3 = chol(inv(Db3));
%     b3 = Db3\db3 + H3'*randn(k,1);
%     %Update Sige
%     nu_3 = nu3 + Nx;
%     d1_3 = d03 + (ys3 - x*b3)'*(ys3 - x*b3);
%     chi_3 = chis_rnd(1,nu_3);
%     s3 = d1_3/chi_3;
%     %Update rho
%     b0_3 = (xpx + s3*iB03)\(x'*(y-treat*double(r>c)) + s3*(iB03*b03));
%     bd_3 = (xpx + s3*iB03)\(x'*(w*(y-treat*double(r>c))) + s3*(iB03*b03));
%     e0_3 = (y - treat*double(r>c)) - x*b0_3;
%     ed_3 = w*(y-treat*double(r>c)) - x*bd_3;
%     epe0_3 = e0_3'*e0_3;
%     eped_3 = ed_3'*ed_3;
%     epe0d_3 = ed_3'*e0_3;
%     p3 = draw_rho(detval,epe0_3,eped_3,epe0d_3,Nx,k,p3,a13,a23);
    %storage of results
    if j == 1
        beta1(j,:) = b1';
    else
        beta1(j,:) = b2(1:end-1,1)';
    end
        sige1(j,1)   = s2;
        rho1(j,1)    = p2;        
        mbeta2(j,:) = mbeta';
        pbeta2(j,:) = pbeta';
        msige2(j,1)   = s2m;
        psige2(j,1)   = s2p;
        bw_left(j,1) = bw_l;
        bw_right(j,1)= bw_r;
%         beta3(j - nomit,:) = [b3' treat];
%         sige3(j - nomit)   = s3;
%         rho3(j - nomit)    = p3;
        treat_save(j,1) = treat;
        
        
        if ( mod(j, noo) ==0 )
            
            tt=j-noo+1:j-1;
            
            subplot(2,3,1),
            t1 = 1:length(wxp_save);
            t2 = 1:length(wxm_save);
            plot(t1,sort(wxp_save(:,1)),'-g',t1,sort(wxp_save(:,2)),'-r');
            xlabel('wxp');
            subplot(2,3,2),
            plot(t2,sort(wxm_save(:,1)),'-g',t2,sort(wxm_save(:,2)),'-r');
            xlabel('wxm');
            subplot(2,3,3),
            plot(tt,rho1(tt,1));
            xlabel('rho1');
            subplot(2,3,4),
            plot(tt,treat_save(tt,1));
            xlabel('treatment  draws');
            subplot(2,3,5),
            plot(tt,beta1(tt,1),tt,beta1(tt,k/2+1));
            xlabel('beta1, theta1');
            subplot(2,3,6),
            plot(tt,pbeta2(tt,1),tt,mbeta2(tt,1));
            xlabel('treatment b1, b2');
            
            drawnow;
        end
%     waitbar(j/ndraw);
end
% close(hwait);

% result.bdraw = beta1;
% result.pdraw = p1;
% result.sdraw = s1;
result.tr = treat_save(nomit+1:end,1);
result.samples = samples;
% result.pbeta = pbeta2;
% result.mbeta = mbeta2;

%         beta1(j,:) = b1';
%         sige1(j,1)   = s1;
%         rho1(j,1)    = p1;
%         mbeta2(j,:) = mbeta';
%         pbeta2(j,:) = pbeta';
%         msige2(j,1)   = s2m;
%         psige2(j,1)   = s2p;
%         bw_left(j,1) = bw_l;
%         bw_right(j,1)= bw_r;
% %         beta3(j - nomit,:) = [b3' treat];
% %         sige3(j - nomit)   = s3;
% %         rho3(j - nomit)    = p3;
%         treat_save(j,1) = treat;


%Calculating partial effects
uiter = 50;
maxorderu = 150;
ny = length(y);
rv = randn(ny,uiter);
tracew = zeros(maxorderu,1);
wjjju = rv;
for jjj = 1:maxorderu
   wjjju = W*wjjju;
   tracew(jjj) = mean(mean(rv.*wjjju));
end
traces = tracew;
traces(1) = 0;
traces(2) = sum(sum(W'.*W))/ny;
trs = [1;traces];
ntrs = length(trs);
trbig = trs';
trbig2 = [trbig(1,2:end) trbig(1,end)];
trmat = [trbig
          trbig2];
bdraws = beta1;
pdraws = rho1;
ree = 0:1:ntrs-1;
rmat = zeros(1,ntrs);
total = zeros(ndraw - nomit,k,ntrs);
direct = zeros(ndraw - nomit,k,ntrs);
indirect = zeros(ndraw - nomit,k,ntrs);

for i = nomit+1:ndraw
    rmat = pdraws(i).^ree;
    for j = 1:k/2
%         beta = [bdraws(i,j) bdraws(i,j+pp)];
        beta = [bdraws(i,j)  bdraws(i,j+k/2)];

        total(i-nomit,j,:) = (beta(1) + beta(2))*rmat;
        direct(i-nomit,j,:) = (beta*trmat).*rmat;
        indirect(i-nomit,j,:) = total(i-nomit,j,:) - direct(i-nomit,j,:);
    end
end

total_treat = zeros(ndraw-nomit,ntrs);
direct_treat = zeros(ndraw-nomit,ntrs);
indirect_treat = zeros(ndraw-nomit,ntrs);

for i = nomit+1:ndraw
    beta = treat_save(i,1);
    total_treat(i-nomit,:) =  beta*rmat;
    direct_treat(i-nomit,:) = (beta*trbig).*rmat;
    indirect_treat(i-nomit,:) = (beta*rmat) - ((beta*trbig).*rmat);
end

total_dist = zeros(ndraw-nomit,k);
direct_dist = zeros(ndraw-nomit,k);
indirect_dist = zeros(ndraw-nomit,k);
for i = 1:k
    total_dist(:,i) = sum(squeeze(total(:,i,:)),2);
    direct_dist(:,i) = sum(squeeze(direct(:,i,:)),2);
    indirect_dist(:,i) = sum(squeeze(indirect(:,i,:)),2);
end



%Setting up Outputs
result.beta = mean(beta1(nomit+1:ndraw,:));
result.rho  = mean(rho1(nomit+1:ndraw,1));
result.sige = mean(sige1(nomit+1:ndraw,1));
result.lower = mbeta2(nomit+1:ndraw,1:2);
result.upper = pbeta2(nomit+1:ndraw,1:2);
result.bdraw = beta1(nomit+1:ndraw,:);
result.pdraw = rho1(nomit+1:ndraw,:);
result.sdraw = sige1(nomit+1:ndraw,:);
result.total_treat_draws = sum(total_treat,2);
result.direct_treat_draws = sum(direct_treat,2);
result.indirect_treat_draws = sum(indirect_treat,2);
result.total_draws = total_dist;  
result.direct_draws = direct_dist;
result.indirect_draws = indirect_dist;
result.total_treat_mean = mean(sum(total_treat,2));
result.direct_treat_mean = mean(sum(direct_treat,2));
result.indirect_treat_mean = mean(sum(indirect_treat,2));
result.total_mean = mean(total_dist);
result.direct_mean = mean(direct_dist);
result.indirect_mean = mean(indirect_dist);
result.samples = samples;
end


function rho = draw_rho(detval,epe0,eped,epe0d,n,k,rho)
% PURPOSE: draws rho-values using griddy Gibbs and inversion
% ---------------------------------------------------
%  USAGE: rho = draw_rho(detval,epe0,eped,epe0d,n,k,rho,a1,a2)
% where info contains the structure variable with inputs 
% and the outputs are either user-inputs or default values
% ---------------------------------------------------
% REFERENCES: LeSage and Pace (2009) Introduction to Spatial Econometrics
% Chapter 5, pp 136-141 on Bayesian spatial regression models.


% written by:
% James P. LeSage, last updated 3/2010
% Dept of Finance & Economics
% Texas State University-San Marcos
% 601 University Drive
% San Marcos, TX 78666
% jlesage@spatial-econometrics.com


nmk = (n-k)/2;
nrho = length(detval(:,1));
iota = ones(nrho,1);

z = epe0*iota - 2*detval(:,1)*epe0d + detval(:,1).*detval(:,1)*eped;
z = -nmk*log(z);
den =  detval(:,2) + z;

% bprior = beta_prior(detval(:,1),a1,a2);
% den = den + log(bprior);
n = length(den);
y = detval(:,1);
adj = max(den);
den = den - adj;
x = exp(den);
% trapezoid rule
isum = sum((y(2:n,1) + y(1:n-1,1)).*(x(2:n,1) - x(1:n-1,1))/2);
z = abs(x/isum);
den = cumsum(z);

rnd = unif_rnd(1,0,1)*sum(z);
ind = find(den <= rnd);
idraw = max(ind);
if (idraw > 0 & idraw < nrho)
rho = detval(idraw,1);
end
end
function [detval,time1] = sar_lndet(ldetflag,W,rmin,rmax,detval,order,iter)
% PURPOSE: compute the log determinant |I_n - rho*W|
% using the user-selected (or default) method
% ---------------------------------------------------
%  USAGE: detval = far_lndet(lflag,W,rmin,rmax)
% where eflag,rmin,rmax,W contains input flags 
% and the outputs are either user-inputs or default values
% ---------------------------------------------------

% written by:
% James P. LeSage, last updated 3/2010
% Dept of Finance & Economics
% Texas State University-San Marcos
% 601 University Drive
% San Marcos, TX 78666
% jlesage@spatial-econometrics.com


% do lndet approximation calculations if needed
if ldetflag == 0 % no approximation
t0 = clock;    
out = lndetfull(W,rmin,rmax);
time1 = etime(clock,t0);
tt=rmin:.001:rmax; % interpolate a finer grid
outi = interp1(out.rho,out.lndet,tt','spline');
detval = [tt' outi];
    
elseif ldetflag == 1 % use Pace and Barry, 1999 MC approximation

t0 = clock;    
out = lndetmc(order,iter,W,rmin,rmax);
time1 = etime(clock,t0);
tt=rmin:.001:rmax; % interpolate a finer grid
outi = interp1(out.rho,out.lndet,tt','spline');
detval = [tt' outi];

elseif ldetflag == 2 % use Pace and Barry, 1998 spline interpolation

t0 = clock;
out = lndetint(W,rmin,rmax);
time1 = etime(clock,t0);
tt=rmin:.001:rmax; % interpolate a finer grid
outi = interp1(out.rho,out.lndet,tt','spline');
detval = [tt' outi];

elseif ldetflag == -1 % the user fed down a detval matrix
    time1 = 0;
        % check to see if this is right
        if detval == 0
            error('sar_g: wrong lndet input argument');
        end
        [n1,n2] = size(detval);
        if n2 ~= 2
            error('sar_g: wrong sized lndet input argument');
        elseif n1 == 1
            error('sar_g: wrong sized lndet input argument');
        end          
end
end
function rn=chis_rnd(nn,v)
% PURPOSE: generates random chi-squared deviates
%---------------------------------------------------
% USAGE:   rchi = chis_rnd(n,v)
% where:   n = a scalar for the size of the vector to be generated
%          v = the degrees of freedom 
% RETURNS: n-vector with mean=v, variance=2*v
% --------------------------------------------------
% SEE ALSO: chis_d, chis_inv, chis_cdf, chis_pdf
% --------------------------------------------------

% This is code by 
%        Anders Holtsberg, 18-11-93
%        Copyright (c) Anders Holtsberg
%   Documentation for the routine compatible was changed
%   to make is compatible with the needs of the
%   Econometrics Toolbox


rn = rgamma(nn,v*0.5)*2;


function x = rgamma(nn,a)
%RGAMMA   Random numbers from the gamma distribution
%
%         x = rgamma(n,a)

% GNU Public Licence Copyright (c) Anders Holtsberg 10-05-2000.

% This consumes about a third of the execution time compared to 
% the Mathworks function GAMRND in a third the number of
% codelines. Yihaaa! (But it does not work with different parameters)
%
% The algorithm is a rejection method. The logarithm of the gamma 
% variable is simulated by dominating it with a double exponential.
% The proof is easy since the log density is convex!
% 
% Reference: There is no reference! Send me an email if you can't
% figure it out.

if any(any(a<=0))
   error('Parameter a is wrong')
end

n = prod(nn);
if length(nn) == 1
   nn(2) = 1;
end

y0 = log(a)-1/sqrt(a);
c = a - exp(y0);
m = ceil(n*(1.7 + 0.6*(min(min(a))<2)));

y = log(rand(m,1)).*sign(rand(m,1)-0.5)/c + log(a);
f = a*y-exp(y) - (a*y0 - exp(y0));
g = c*(abs((y0-log(a))) - abs(y-log(a)));
reject = (log(rand(m,1)) + g) > f;
y(reject) = [];
% x = zeros(n,1);
if length(y) >= n
   x = exp(y(1:n));
else
   tmp = rgamma(n - length(y), a);
   x = [exp(y) 
        tmp];
end
end
end
function [wyp, wym, wxp, wxm, dplus, dminus] = weights(y,r,c,bw_l,bw_r)
    rc = r-c;
    d  = rc>0;
    rp = rc(rc>=0,1);
    rm = rc(rc<0,1);
    yp = y(rc>=0,1);
    ym = y(rc<0,1);
    dp = d(rc>=0,1);
    dm = d(rc<0,1);
    %on plus side
    weight = (rp<=bw_r).*(1-rp/bw_r);
    sweight = sqrt(weight((rp<=bw_r),1));
    wxp = [sweight,rp(abs(rp)<=bw_r,:).*(sweight*ones(1,1))];
    wyp = yp(rp<=bw_r).*sweight;
    dplus = dp(rp<=bw_r).*sweight;
    %on minus side
    weight = ((abs(rm)<=bw_l)).*(1-abs(rm)/bw_l);
    sweight = sqrt(weight(((abs(rm)<=bw_l)),1));
    wxm = [sweight,rm(abs(rm)<=bw_l,:).*(sweight*ones(1,1))];
    wym = ym((abs(rm)<=bw_l)).*sweight;
    dminus = dm(abs(rm)<=bw_l).*sweight;
end
function output = rddmmse(X,Y,c)

%--------------------------------------------------------------------------
% rddmmse.m: compute the sharp regression discontinuity estimator
%--------------------------------------------------------------------------
%
% DESCRIPTION:
%   Estimate the sharp regression discontinuity estimator with the
%   bandwidths proposed by Arai and Ichimura (2015).
%
% REFERENCE:
%   Arai and Ichimura (2015). "Simultaneous Selection of Optimal Bandwidths 
%   for the Sharp Regression Discontinuity Estimator," CIRJE Discussion Paper,
%   University of Tokyo, CIRJE-F-984.
%
% 
% INPUT ARGUMENTS
%   X: Imput vector (n by 1) of the forcing variable 
%   Y: Vector of the outcome variable (n by 1) 
%   c: Evaluation (discontinuity) point (scalar)
%
% OUTPUT ARGUMENTS:
%   tau:    sharp RD estimate
%   se:     standard error
%   h:      bandwidths, the first element for the right of the cut-off
%                   the second element for the left of the cut-off
%   tval:   bias-corrected t-value
%      ci:  bias-corrected confidence interval, See Section 4.4 of
%               Fan and Gijbels (1996)
%       n:  the number of observations with nonzero weight,
%               the first element for the right of the cut-off,
%               the second element for the left of the cut-off
% 
% AUTHOR:
%   This code was written by Yoichi Arai.
%   I will appreciate any feedbacks.
%
% Note:
%   Optimization toolbox is necessary.
% 
%
% 
% DEVELOPMENT: 
%    3 April 2014: Original version of rddmmse.m written.
%    4 September 2015: Minor changes in Step 2 and 3 (See Arai and Ichimura
%    (2015)
%--------------------------------------------------------------------------

% Check input arguments
if nargin < 3
    error('more input arguments needed.');
end

if nargin > 3
    error('Too many arguments.');
end

% Step 1: Estimation of density and its first derivative
 % Density Estimation
 n = size(X,1);
 lambda = sqrt(var(X));
 a_pilot = 2.34*lambda/n^(1/5);
 f_k = kdensity_est(X,c,a_pilot); % kernel = 2, Epanechnikov
 h1_pilot = lambda*(105*16*sqrt(pi)/15/n)^(1/7);
 f1_pilot = density1d_est(X,c,h1_pilot);
 
% Step 2: Get pilot bandwidths to estimate derivatives
 % 4th order Polynomial regression
 Xbin_lt_x = (X<c);           % binary indicator for observations less than x
 Xbin_bt_x = (X>=c);          % bigger than x
 num_X_bt_x = sum(Xbin_bt_x); % # of observations of X bigger than x, N+
 num_X_lt_x = sum(Xbin_lt_x); % # of observations of X less than x, N- in p.16
 Xx = X-c;
 Xx2 = Xx.^2;
 Xx3 = Xx.*Xx2;
 Xx4 = Xx2.^2;
 nones = ones(n,1);
 data_mat = [Y nones Xx Xx2 Xx3 Xx4];
 data_mat_p = data_mat(Xbin_bt_x==1,:); % data_mat only for positive X
 data_mat_m = data_mat(Xbin_lt_x==1,:); % data_mat only for negative X

 YT_pc = data_mat_p(:,1);
 TT_pc = data_mat_p(:,2:6);
 lambda_pc_4_4g = TT_pc\YT_pc;
 m41_4 = 24*lambda_pc_4_4g(5,1);
 
 YT_mc = data_mat_m(:,1);
 TT_mc = data_mat_m(:,2:6);
 lambda_mc_4_4g = TT_mc\YT_mc;
 m42_4 = 24*lambda_mc_4_4g(5,1);

 vy_p_plugin_4g = sum((YT_pc-TT_pc*lambda_pc_4_4g).^2)/(num_X_bt_x-5);
 vy_m_plugin_4g = sum((YT_mc-TT_mc*lambda_mc_4_4g).^2)/(num_X_lt_x-5);

 % Pilot bandwidths for 3rd order LPR to estimate 3rd derivative with uniform kernel
 C33K_U = Const_LPR_U_est(3,3);
 h3_p = C33K_U*(vy_p_plugin_4g/f_k/m41_4^2)^(1/9)/num_X_bt_x^(1/9);
 h3_m = C33K_U*(vy_m_plugin_4g/f_k/m42_4^2)^(1/9)/num_X_lt_x^(1/9);
 % Pilot bandwidths for 3rd order LPR to estimate 2nd derivative with uniform kernel
 C23K_U = Const_LPR_U_est(2,3);
 h23_p = C23K_U*(vy_p_plugin_4g/f_k/m41_4^2/num_X_bt_x)^(1/9); % possible to use either v_pilot_24 or v_pilot_23
 h23_m = C23K_U*(vy_m_plugin_4g/f_k/m42_4^2/num_X_lt_x)^(1/9);
 

% Step 3: Estimate the second and the third derivatives
 % 3rd order LPR for the third derivatives    
 Xbin_bet_xxh = Xbin_bt_x.*(X<=c+h3_p);     % index for X in [x, x+h3_p]
 data_mat_4_p = data_mat(Xbin_bet_xxh==1,:);
 YT_p = data_mat_4_p(:,1);
 T_p = data_mat_4_p(:,2:5);
 lambda_p_all_33_4 = T_p\YT_p;

 Xbin_bet_xhx = (X>=c-h3_m).*Xbin_lt_x;     % index for X in [x-h3_m, x]
 data_mat_4_m = data_mat(Xbin_bet_xhx==1,:);
 YT_m = data_mat_4_m(:,1);
 T_m = data_mat_4_m(:,2:5);
 lambda_m_all_33_4 = T_m\YT_m;   

 m3c_p = 6*lambda_p_all_33_4(4,1);
 m3c_m = 6*lambda_m_all_33_4(4,1);

 % 3rd order LPR for the second derivatives
 Xbin_bet_xxh = Xbin_bt_x.*(X<=c+h23_p);     % index for X in [x, x+h3_p]
 data_mat_4_p = data_mat(Xbin_bet_xxh==1,:);
 YT_p = data_mat_4_p(:,1);
 T_p = data_mat_4_p(:,2:5);
 lambda_p_all_23 = T_p\YT_p;
 m2c_p_2 = 2*lambda_p_all_23(3,1);

 YT_p_hat = T_p*lambda_p_all_23;
 vy_p_4_2 = sum((YT_p - YT_p_hat).^2)/(size(YT_p,1)-4);

 Xbin_bet_xhx = (X>=c-h23_m).*Xbin_lt_x;     % index for X in [x-h3_m, x]
 data_mat_4_m = data_mat(Xbin_bet_xhx==1,:);
 YT_m = data_mat_4_m(:,1);
 T_m = data_mat_4_m(:,2:5);
 lambda_m_all_23 = T_m\YT_m;
 m2c_m_2 = 2*lambda_m_all_23(3,1);

 YT_m_hat = T_m*lambda_m_all_23;
 vy_m_4_2 = sum((YT_m - YT_m_hat).^2)/(size(YT_m,1)-4);

 m2_pilot = [m2c_p_2 m2c_m_2];
 m3_pilot = [m3c_p m3c_m];
 v_pilot = [vy_p_4_2 vy_m_4_2];

% Step 4: Get bandwidths that minimize MMSE
 options = optimset('Algorithm','interior-point','Display','off');
 Qn_value_2 = @(theta)Qn(theta,f_k,f1_pilot,m2_pilot,m3_pilot,v_pilot,n);

 % Lower bound for AI search area
 NN_min_AI = 3;
 X_lt_x = data_mat_m(:,3);
 X_bt_x = data_mat_p(:,3);
 sorted_X_bt_x = sort(X_bt_x);               % ascending order
 sorted_X_lt_x_des_abs = abs(sort(X_lt_x,'descend'));
 NN_p_AI = sorted_X_bt_x(NN_min_AI,1);
 NN_m_AI = sorted_X_lt_x_des_abs(NN_min_AI,1);
 h_lb = [NN_p_AI,NN_m_AI];
 % Upper bound for AI search area
 hu1 = sorted_X_bt_x(num_X_bt_x,1);
 hu2 = sorted_X_lt_x_des_abs(num_X_lt_x,1);
 h_ub = [hu1,hu2];

 h_mat_local = zeros(10,2);
 min_val = zeros(10,1);
 h1_range = sorted_X_bt_x(num_X_bt_x,1) - sorted_X_bt_x(3,1);
 h2_range = sorted_X_lt_x_des_abs(num_X_lt_x,1) - sorted_X_lt_x_des_abs(3,1);
 
 
 for k=1:10
 for j=1:10
     h_ini_1 = sorted_X_bt_x(3,1) + 0.01*h1_range*j;
     h_ini_2 = sorted_X_lt_x_des_abs(3,1) + 0.1*h2_range*k;
     h_ini = [h_ini_1 h_ini_2];
     [theta_2, fval_2] = fmincon(Qn_value_2,h_ini,[],[],[],[],h_lb,h_ub,[],options);
     h_mat_local(10*(k-1)+j,:) = theta_2;
     min_val(10*(k-1)+j,1) = fval_2;
     
 end
 end
 [~, min_index] = min(min_val);
 h_opt = h_mat_local(min_index,:);
 [tau,se,~,~,n_p,n_m] = srdd_est(Y,X,c,h_opt(1),h_opt(2));
 
% Compute bias correction terms b11 and b12
Xbin_bt_xh_opt = (X >= c-h_opt(2));
Xbin_lt_xh_opt = (X <= c+h_opt(1));
Xbin_bet_xxh_opt = Xbin_bt_x.*Xbin_lt_xh_opt;    % index for X in [x-h,x]
Xbin_bet_xhx_opt = Xbin_bt_xh_opt.*Xbin_lt_x;    % index for X in [x,x+h]

data_trim_mat_p_opt = data_mat(Xbin_bet_xxh_opt==1,:);
T_p_opt = data_trim_mat_p_opt(:,2:3);
Kvec_p_opt = (1-abs(data_trim_mat_p_opt(:,3))/h_opt(1)).*(abs(data_trim_mat_p_opt(:,3))<h_opt(1));
W_p_opt = diag(Kvec_p_opt);

data_trim_mat_m_opt = data_mat(Xbin_bet_xhx_opt==1,:);
T_m_opt = data_trim_mat_m_opt(:,2:3);
Kvec_m_opt = (1-abs(data_trim_mat_m_opt(:,3))/h_opt(2)).*(abs(data_trim_mat_m_opt(:,3))<h_opt(2));
W_m_opt = diag(Kvec_m_opt);
 
Sn01_p_tilde_inv_opt = invpd(T_p_opt'*W_p_opt*T_p_opt);
cn2_p_tilde_opt = data_trim_mat_p_opt(:,4:5)'*Kvec_p_opt;
cn3_p_tilde_opt = data_trim_mat_p_opt(:,5:6)'*Kvec_p_opt;
b11_vec_opt = m2c_p_2*Sn01_p_tilde_inv_opt*cn2_p_tilde_opt/2 ...
    + m3c_p*Sn01_p_tilde_inv_opt*cn3_p_tilde_opt/6;

Sn01_m_tilde_inv_opt = invpd(T_m_opt'*W_m_opt*T_m_opt);
cn2_m_tilde_opt = data_trim_mat_m_opt(:,4:5)'*Kvec_m_opt;
cn3_m_tilde_opt = data_trim_mat_m_opt(:,5:6)'*Kvec_m_opt;
b12_vec_opt = m2c_m_2*Sn01_m_tilde_inv_opt*cn2_m_tilde_opt/2 ...
    + m3c_m*Sn01_m_tilde_inv_opt*cn3_m_tilde_opt/6;

bias_opt = b11_vec_opt(1,1)-b12_vec_opt(1,1);

% Construct bias-corrected t-value
tval = (tau-bias_opt)/se;

output.tau = tau;
output.se =se;
output.h = h_opt;
output.tval = tval;
output.ci = [tau-bias_opt-1.96*se tau-bias_opt+1.96*se];
output.n = [n_p n_m];
%------------------------------------------------------------------
% End function rddmmse
%------------------------------------------------------------------
end
function fhat = kdensity_est(X,c,h)

% Purpose: estimate density at a point with Epanechnikov kernel
% --------------------------------------------
% Usage: results = density_est
%       X = a column vector of scalar random variables
%       c = point to estimate
%       h = bandwidth
% --------------------------------------------

nobs = size(X,1);
%fhat = 0;
A_h = (X-c)/h;
ExpA_h = 3*sum((abs(A_h)<1).*(1-A_h.^2))/4;
fhat = ExpA_h/nobs/h;

%------------------------------------------------------------------
% End private function density_est
%------------------------------------------------------------------
end
function f1hat = density1d_est(X,c,h)

% Purpose: estimate 1st derivative density at a point
% --------------------------------------------
% Usage: results = density_est
%       X = a column vector of scalar random variables
%       c = point to estimate
%       h = bandwidth
% --------------------------------------------

nobs = size(X,1);
A_h = (c-X)/h;
B_h = A_h.^2;
k2 =1/5;
const = 3/4;
ExpA_h = (-1)*const*sum(A_h.*(1-B_h).*(abs(A_h)<1))/k2;
f1hat = ExpA_h/nobs/h^2;

%------------------------------------------------------------------
% End private function density1d_est
%------------------------------------------------------------------
end
function [C_vp] = Const_LPR_U_est(v,p)

% S = (\mu_{j+k}), 0<=j,k<=p, defined in FG(96,p61), p=6
S5 = [1/2 1/4 1/6 1/8 1/10;
       1/4 1/6 1/8 1/10 1/12;
       1/6 1/8 1/10 1/12 1/14;
       1/8 1/10 1/12 1/14 1/16;
       1/10 1/12 1/14 1/16 1/18];
% Sstar = (\nu_{j+k}), 0<=j,k<=p, defined in FG(96,p61), p=6
K5 = [1/4 1/8 1/12 1/16 1/20;
        1/8 1/12 1/16 1/20 1/24;
        1/12 1/16 1/20 1/24 1/28;
        1/16 1/20 1/24 1/28 1/32;
        1/20 1/24 1/28 1/32 1/36];
        
% Vector of \mu_j,
% [\mu_1, \mu_2, ..., \mu_{13}]
mu_vec = [1/4 1/6 1/8 1/10 1/12 1/14 1/16 1/18 1/20];

Sinv_p = invpd(S5(1:p+1,1:p+1));

% int t^{p+1}K_{v,p}^*(t)dt
int_tp_K = Sinv_p(v+1,:)*mu_vec(1,(p+1):(2*p+1))';

% \int_K_{v,p}^{*2}(t)dt
int_K2 = Sinv_p(v+1,:)*K5(1:p+1,1:p+1)*Sinv_p(v+1,:)';

C_vp = ((factorial(p+1)^2*(2*v+1)*int_K2)/(2*(p+1-v)*int_tp_K^2))^(1/(2*p+3));

%------------------------------------------------------------------
% End private function Const_LPR_U_est
%------------------------------------------------------------------
end
function [Qn] = Qn(h,f_pilot,f1_pilot,m2_pilot,m3_pilot,v_pilot,n)

% matrix for constants depending on kernel
% [b1 c1 c2 v]
k_const = [-0.1 -0.1 -0.08 24/5];

b2 = zeros(1,2);
m2foverf12 = zeros(1,2);
for j=1:2
    m2foverf12(j) = m2_pilot(j)*f1_pilot/f_pilot/2;
    b2(j) = (-1)^(j+1)*(k_const(2)*(m2foverf12(j) + m3_pilot(j)/6)...   
         - k_const(3)*m2foverf12(j));
end
    
Qn = k_const(1)^2*(m2_pilot(1)*h(1)^2 - m2_pilot(2)*h(2)^2)^2/4 ...
+(b2(1)*h(1)^3 - b2(2)*h(2)^3)^2 ...
+ k_const(4)*(v_pilot(1)/h(1) + v_pilot(2)/h(2))/f_pilot/n;

%------------------------------------------------------------------
% End private function Qn
%------------------------------------------------------------------
end

function [h,flag] = ikbw (Y,X,c)
%estimates optimal bandwidth for sharp RD with triangular kernel following
%Imbens and Kalyanaraman(2012) - IK
%
%
%Marinho Bertanha, Wei Qian
%This version: 11/07/2018
%optimal bandwidth algorithm based on rd_optbandwidth.m by IK
%
%
%
%
%Input
%
%Y - (n x 1) sample of Y, outcome variable
%X - (n x 1) sample of X, forcing variable
%c - scalar cutoff


%
%Output
%h - optimal bandwidth for sharp RD according to IK
%
%flag (integer) equals to the number of matrices that were not invertible;
%in this case, a pseudo-inverse is used; 0 o.w.;
%


flag=0;
nmin=5;

n=length(Y);
nl=sum(X<c);
nr=sum(X>=c);




%**************************************************************************
%density and side limits of variance

hx=1.84*std(X)*(n^(-1/5)); %Silverman's bandwidth for uniform kernel and Normal X



%restricts sample to usable observations
Ylh=Y((X<c)&(c-hx<X));
Yrh=Y((X>=c)&(X<c+hx));
 


nlh=sum((X<c)&(c-hx<X));
nrh=sum((X>=c)&(X<c+hx));

if min(nlh,nrh)<nmin
    display('1: Not enough observations around the cutoff! Try a larger bandwidth or sample.');
    flag=flag+1;
    return;
end

%density
fx=(nlh+nrh)/(2*n*hx);

%variance
syl=(1/(nlh-1))*sum((Ylh-mean(Ylh)).^2);
syr=(1/(nrh-1))*sum((Yrh-mean(Yrh)).^2);



%**************************************************************************
%second derivatives of m_y
XX = [ones(n,1) (X>=c) (X-c) (X-c).^2 (X-c).^3];

if rank(XX'*XX)<size(XX'*XX,1)
    flag=flag+1;
    my3=[0 0 0 0 6]*pinv(XX'*XX)*(XX'*Y);

else
    my3=[0 0 0 0 6]*(eye(5)/(XX'*XX))*(XX'*Y);

end



hy2l = 3.56*((syl/(fx*my3^2))^(1/7))*(nl^(-1/7));
hy2r = 3.56*((syr/(fx*my3^2))^(1/7))*(nr^(-1/7));



Ylh=Y((X<c)&(c-hy2l<X));
Yrh=Y((X>=c)&(X<c+hy2r));

Xylh=X((X<c)&(c-hy2l<X));
Xyrh=X((X>=c)&(X<c+hy2r));

nylh=sum((X<c)&(c-hy2l<X));
nyrh=sum((X>=c)&(X<c+hy2r));

if min(nylh,nyrh)<nmin
    display('2: Not enough observations around the cutoff! Try a larger bandwidth or sample.');
    flag=flag+1; 
    return;
end

XXylh = [ones(nylh,1) (Xylh-c) (Xylh-c).^2];

if rank(XXylh'*XXylh)<size(XXylh'*XXylh,1)
    flag=flag+1;
    my2l=[0 0 2]*pinv(XXylh'*XXylh)*(XXylh'*Ylh);

else
    my2l=[0 0 2]*(eye(3)/(XXylh'*XXylh))*(XXylh'*Ylh);    

end


XXyrh = [ones(nyrh,1) (Xyrh-c) (Xyrh-c).^2];

if rank(XXyrh'*XXyrh)<size(XXyrh'*XXyrh,1)
    flag=flag+1;
    my2r=[0 0 2]*pinv(XXyrh'*XXyrh)*(XXyrh'*Yrh);

else
    my2r=[0 0 2]*(eye(3)/(XXyrh'*XXyrh))*(XXyrh'*Yrh);

end


%**************************************************************************
%regularization terms
rl = 2160*syl/(nylh*(hy2l^4));
rr = 2160*syr/(nyrh*(hy2r^4));

%**************************************************************************


h = 3.4375*(((syl+syr)/(fx*((my2r-my2l)^2+(rl+rr))))^(1/5))*(n^(-1/5));




end


function [sigma] = match(W,X)
%#codegen
m     = length(W);
y_ave = zeros(m,1);
for k=1:m
   diffx      = abs(X - X(k,:)); 
   m_group    = sort(diffx);
   m_group    = m_group([true;diff(m_group(:))>0]);
   m_group    = m_group(2:4,:);
   ind        = (diffx == m_group(1,1) | diffx == m_group(2,1) | diffx == m_group(3,1));
   y_ave(k,1) = mean(W(ind));
end
sigma = .75*(W - y_ave).^2;

end

function [h,b]=rdbwselect_CCT(w,x,threshold)

% matlab replication of rdbwselect.R/version 0.5 06Jun2014

% INPUTS
% w is the outcome
% x is the scalar forcing variable
% threshold is the value of the threshold for the forcing variable

% OUTPUT	
% h_opt is bandwidth


kerf        = @(z)(abs(z) <= 1).*(1 - abs(z));
c           = threshold;
N           = length(w);% number of observations
x           = x-c;% forcing variable relative to threshold
obs_left    = (x < 0);%Observations to left of threshold in each sample
X_l         = x(obs_left,:);
W_l         = w(obs_left,:);
n           = sum(obs_left);%Number of observation to left of threshold in each sample
obs_right   = (x >= 0);%Observations to right of threshold in each sample
X_r         = x(obs_right,:);
W_r         = w(obs_right,:);
m           = sum(obs_right);

%----STEP 0

%------Calculating v_n

S_x              = std(x);
IQR_x            = iqr(x);
temp             = [S_x (IQR_x/1.349)];
v_n              = 2.58*min(temp)*N^(-1/5);

%------Calculating c_n
%------Calculating V_33_plus(v_n) + V_33_minus(v_n)

X_plus_3         = [ones(m,1) (X_r) (X_r).^2 (X_r).^3];
W_plus           = diag(kerf(X_r./v_n)./v_n);
Psi_plus_3       = X_plus_3'*W_plus*X_plus_3;
temp             = (X_r < v_n);
sigma_plus       = zeros(m,1);
sigma_plus(temp,:)  = match(W_r(temp,:),X_r(temp,:));
sigma_plus       = diag(sigma_plus);
Phi_plus_33      = (X_plus_3'*W_plus*sigma_plus*W_plus*X_plus_3);
V_plus_33        = ([0 0 0 1]*inv(Psi_plus_3)*Phi_plus_33*inv(Psi_plus_3)*([0 0 0 1].'));

X_minus_3        = [ones(n,1) (X_l) (X_l).^2 (X_l).^3];
W_minus          = diag(kerf(X_l./v_n)./v_n);
Psi_minus_3      = X_minus_3'*W_minus*X_minus_3;
temp             = (X_l > -v_n);
sigma_minus      = zeros(n,1);
sigma_minus(temp,:)  = match(W_l(temp,:),X_l(temp,:));
sigma_minus      = diag(sigma_minus);
Phi_minus_33     = (X_minus_3'*W_minus*sigma_minus*W_minus*X_minus_3);
V_minus_33       = ([0 0 0 1]*inv(Psi_minus_3)*Phi_minus_33*inv(Psi_minus_3)*([0 0 0 1].'));

V_33             = V_plus_33 + V_minus_33;

%------Calculating B_33

B_33             = 1.777778;

%------Calculating global estimates gamma_4_plus + gamma_4_minus

X_plus_4         = [ones(m,1) X_r (X_r).^2 (X_r).^3 (X_r).^4];  
gamma_plus_4     = inv(X_plus_4'*X_plus_4)*X_plus_4'*W_r;

X_minus_4        = [ones(n,1) X_l (X_l).^2 (X_l).^3 (X_l).^4];  
gamma_minus_4    = inv(X_minus_4'*X_minus_4)*X_minus_4'*W_l ;

%------Calculating c_n using above

Num    = 7*N*v_n^(7)*V_33;
Den    = 2*B_33^2*(gamma_plus_4(5,1) - (-1)^(2)*gamma_minus_4(5,1))^2;
C_012  = Num/Den;
c_n    = C_012^(1/9)*N^(-1/9);

%----STEP 1

%------Calculating b_n
%------Calculating V_22_plus(v_n) + V_22_minus(v_n)

X_plus_2         = [ones(m,1) X_r X_r.^2];
W_plus           = diag(kerf(X_r./v_n)./v_n);
Psi_plus_2       = X_plus_2'*W_plus*X_plus_2;
temp             = (X_r < c_n);
sigma_plus       = zeros(m,1);
sigma_plus(temp,:)  = match(W_r(temp,:),X_r(temp,:));
sigma_plus       = diag(sigma_plus);
Phi_plus_22      = (X_plus_2'*W_plus*sigma_plus*W_plus*X_plus_2);
V_plus_22        = 2^2*([0 0 1]*inv(Psi_plus_2)*Phi_plus_22*inv(Psi_plus_2)*([0 0 1].'))./(N*v_n^(4));
V_plus_22        = V_plus_22*(N*v_n^(4))/4;

X_minus_2        = [ones(n,1) X_l X_l.^2];
W_minus          = diag(kerf(X_l./v_n)./v_n);
Psi_minus_2      = X_minus_2'*W_minus*X_minus_2;
temp             = (X_l > -c_n);
sigma_minus      = zeros(n,1);
sigma_minus(temp,:)  = match(W_l(temp,:),X_l(temp,:));
sigma_minus      = diag(sigma_minus);
Phi_minus_22     = (X_minus_2'*W_minus*sigma_minus*W_minus*X_minus_2);
V_minus_22       = 2^2*([0 0 1]*inv(Psi_minus_2)*Phi_minus_22*inv(Psi_minus_2)*([0 0 1].'))./(N*v_n^(4));
V_minus_22       = V_minus_22*(N*v_n^(4))/4;

V_22             = V_plus_22 + V_minus_22;

%------Calculating B_22

B_22             = 1.285714;

%------Calculating V_33_plus(c_n) + V_33_minus(c_n)

X_plus_3         = [ones(m,1) X_r X_r.^2 X_r.^3];  
W_plus           = diag(kerf(X_r./c_n)./c_n);
Psi_plus_3       = X_plus_3'*W_plus*X_plus_3;
Phi_plus_33      = (X_plus_3'*W_plus*sigma_plus*W_plus*X_plus_3);
V_plus_33        = 6^2*([0 0 0 1]*inv(Psi_plus_3)*Phi_plus_33*inv(Psi_plus_3)*([0 0 0 1].'))./(N*c_n^(6));
V_plus_33        = V_plus_33*(N*c_n^(6))/36;

X_minus_3        = [ones(n,1) X_l X_l.^2 X_l.^3]; 
W_minus          = diag(kerf(X_l./c_n)./c_n);
Psi_minus_3      = X_minus_3'*W_minus*X_minus_3;
Phi_minus_33     = (X_minus_3'*W_minus*sigma_minus*W_minus*X_minus_3);
V_minus_33       = 6^2*([0 0 0 1]*inv(Psi_minus_3)*Phi_minus_33*inv(Psi_minus_3)*([0 0 0 1].'))./(N*c_n^(6));
V_minus_33       = V_minus_33*(N*c_n^(6))/36;

V_33_c           = V_plus_33 + V_minus_33;

%------Calculating estimates beta_3_plus(c_n) + beta_3_minus(c_n)
  
X_plus_3         = [ones(m,1) X_r X_r.^2 X_r.^3];
W_plus           = diag(kerf(X_r./c_n)./c_n);
Psi_plus_3       = X_plus_3'*W_plus*X_plus_3;
beta_plus_3      = (inv(Psi_plus_3)*X_plus_3'*W_plus*W_r);
 
X_minus_3        = [ones(n,1) X_l X_l.^2 X_l.^3];
W_minus          = diag(kerf(X_l./c_n)./c_n);
Psi_minus_3      = X_minus_3'*W_minus*X_minus_3;
beta_minus_3     = (inv(Psi_minus_3)*X_minus_3'*W_minus*W_l);

%------Calculating b_n using above

Num     = 5*N*v_n^(5)*V_22;
Den     = 2*B_22^2*(([0 0 0 1]*beta_plus_3 - (-1)^3*[0 0 0 1]*beta_minus_3)^2 + 3*V_33_c);
C_012_2 = Num/Den;
b_n     = C_012_2^(1/7)*N^(-1/7);

%----STEP 2

%------Calculating V_01_plus(v_n) + V_01_minus(v_n)

X_plus_1       = [ones(m,1) X_r];
W_plus         = diag(kerf(X_r./v_n)./v_n);
Psi_plus_1     = X_plus_1'*W_plus*X_plus_1;
temp             = (X_r < v_n);
sigma_plus       = zeros(m,1);
sigma_plus(temp,:)  = match(W_r(temp,:),X_r(temp,:));
sigma_plus       = diag(sigma_plus);
Phi_plus_1     = (X_plus_1'*W_plus*sigma_plus*W_plus*X_plus_1);
V_plus_01      = ([1 0]*inv(Psi_plus_1)*Phi_plus_1*inv(Psi_plus_1)*[1 0]')/N;
V_plus_01      = V_plus_01*N;

X_minus_1      = [ones(n,1) X_l];
W_minus        = diag(kerf(X_l./v_n)./v_n);
Psi_minus_1    = X_minus_1'*W_minus*X_minus_1;
temp             = (X_l > -v_n);
sigma_minus      = zeros(n,1);
sigma_minus(temp,:)  = match(W_l(temp,:),X_l(temp,:));
sigma_minus      = diag(sigma_minus);
Phi_minus_1    = (X_minus_1'*W_minus*sigma_minus*W_minus*X_minus_1);
V_minus_01     = ([1 0]*inv(Psi_minus_1)*Phi_minus_1*inv(Psi_minus_1)*[1 0]')/N;
V_minus_01     = V_minus_01*N;

V_01           = V_plus_01 + V_minus_01 ;

%------Calculating B_01

B_01           = -0.1;

%------Calculating estimates beta_2_plus(b_n) + beta_2_plus(b_n)

X_plus_2       = [ones(m,1) X_r X_r.^2];
W_plus         = diag(kerf(X_r./b_n)./b_n);
Psi_plus_2     = X_plus_2'*W_plus*X_plus_2;
beta_plus_2    = (inv(Psi_plus_2)*X_plus_2'*W_plus*W_r);

X_minus_2      = [ones(n,1) X_l X_l.^2];
W_minus        = diag(kerf(X_l./b_n)./b_n);
Psi_minus_2    = X_minus_2'*W_minus*X_minus_2;
beta_minus_2   = (inv(Psi_minus_2)*X_minus_2'*W_minus*W_l);

%------Calculating V_22_plus(b_n) + V_22_minus(b_n)

X_plus_2       = [ones(m,1) X_r X_r.^2];
W_plus         = diag(kerf(X_r./b_n)./b_n);
temp             = (X_r < b_n);
sigma_plus       = zeros(m,1);
sigma_plus(temp,:)  = match(W_r(temp,:),X_r(temp,:));
sigma_plus       = diag(sigma_plus);
Psi_plus_2     = X_plus_2'*W_plus*X_plus_2;
Phi_plus_22    = (X_plus_2'*W_plus*sigma_plus*W_plus*X_plus_2);
V_plus_22      = ([0 0 1]*inv(Psi_plus_2)*Phi_plus_22*inv(Psi_plus_2)*[0 0 1]');

X_minus_2      = [ones(n,1) X_l X_l.^2]; 
W_minus        = diag(kerf(X_l./b_n)./b_n);
Psi_minus_2    = X_minus_2'*W_minus*X_minus_2;
temp             = (X_l > -b_n);
sigma_minus      = zeros(n,1);
sigma_minus(temp,:)  = match(W_l(temp,:),X_l(temp,:));
sigma_minus      = diag(sigma_minus);
Phi_minus_22   = (X_minus_2'*W_minus*sigma_minus*W_minus*X_minus_2);
V_minus_22     = ([0 0 1]*inv(Psi_minus_2)*Phi_minus_22*inv(Psi_minus_2)*[0 0 1]');

V_22_b         = V_plus_22 + V_minus_22;

%------Calculating h_n using above

Num   = N*v_n*V_01;
Den   = 4*B_01^2*(([0 0 1]*beta_plus_2 - (-1)^(2)*[0 0 1]*beta_minus_2)^2 + 3*V_22_b);
C_01  = Num/Den;
h_opt = C_01^(1/5)*N^(-1/5);
h     = h_opt;
b     = b_n;


end

